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Abstract. We use superoperator representation of quantum kinetic equation to 
develop nonequilibrium perturbation theory for inelastic electron current through a 
quantum dot. We derive Lindblad type kinetic equation for an embedded quantum 
dot (i.e. a quantum dot connected to Lindblad dissipators through a buffer zone). The 
kinetic equation is converted to non-Hermitian field theory in Liouville-Fock space. 
The general nonequilibrium many-body perturbation theory is developed and applied 
to the quantum dot with electron-vibronic and electron-electron interactions. Our 
perturbation theory becomes equivalent to Keldysh nonequilibrium Green's functions 
perturbative treatment provided that the buffer zone is large enough to alleviate the 
problems associated with approximations of the Lindblad kinetic equation. 



1. Introduction 

Study of the electron transport through nanoscopic systems remains one of the most 
active areas of contemporary condensed matter physics. Most of the theoretical 
research has been done so far with the use of Keldysh nonequilibrium Green's functions 
(NEGF) pQ and scattering theory based approaches [2] . NEGF applications to electron 
transport were pioneered by Caroli et al.[3] in early 1970s. Keldysh NEGF become 
particularly useful in the development of systematic perturbation theories for electron- 
vibronic and electron-electron interactions in the current-carrying nanosystem. In 
particular, nonequilibrium effects originated from electron-vibration coupling have 
attracted a lot of attention recently because of their importance in single-molecule 
electronics jU EJ EJ El [8] . Various kinds of perturbation theories to deal with electronic 
correlations have been also recently developed (9J EH dU [T21 LTHJ [Hj. 

The electron transport through the system of interacting electrons (either with 
themselves or with some vibrational fields) involves two different energy scales: 
One energy scale is related to the tunneling coupling between the nanosystem and 
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macroscopic leads and the second one is the strength of the interactions inside the 
nanosystems. NEGF usually treats the tunneling interaction exactly, but it has to rely 
on various types of perturbative calculations to account for correlations. On the other 
hand, the approaches based on kinetic equations are able to treat the correlations inside 
the nanosystem very accurately (even exactly in the case of simple model systems) 
but the tunneling part is usually taken into account in the second or sometimes higher 
orders perturbation theory [151 HS1 El HH EEH1 EH 121]. This immediately rules out the 
application of kinetic equations to the one of the most interesting transport regimes 
when there is no energy scale separation between coupling to the electrode and the 
correlations in the the systems (in other words, to the case when the tunneling time 
for electron becomes comparable with the characteristic time for the development of 
correlations in the dot). 

Our approach to the use of kinetic equations for electron transport is different and 
will be elaborated in details in the Sec. [2j We begin with relatively simple kinetic 
equation of the Lindblad type, but we make it exact for the nonequilibrium steady state 
by the introduction of the finite buffer zones between the quantum dot and macroscopic 
leads (so called embedding of the quantum dot) [221 E31 El]- To fully link transport 
kinetic equations with the many-body methods we transform it to Liouville-Fock (or 
super-Fock) space and it becomes equivalent to effective non-Hermitian field theory 
with the right vacuum vector, which corresponds to nonequilibrium steady state density 
matrix. This combination of the embedding and the use of Liouville-Fock space enables 
us to overcome the usual limitations of the kinetic equation based approaches. The 
main goal of the paper is mostly methodological. Namely, we develop nonequilibrium 
perturbation theory in terms of electron-vibronic and electron-electron interaction and 
test our theory against the NEGF results obtained for out of equilibrium local Holstein 
and Anderson models. 

The rest of the paper is organized as follows. In Sec. El we derive the Lindblad 
equation for embedded quantum dot and discuss the underlying approximations. In 
Sec. El we also describe superoperator formalism and convert the kinetic equation to non- 
Hermitian field theory in Liouville-Fock space. Section [3] presents the main equations 
of nonequilibrium many-body perturbation theory, applications to local Holstein and 
Anderson models, and comparison with NEGF. Conclusions are given in Sec. HI We 
use natural units throughout the paper: h = k B = \e\ = 1, where — \e\ is the electron 
charge. 

2. Lindblad kinetic equation for embedded quantum system in 
Liouville-Fock space 

2.1. Lindblad kinetic equation for embedded quantum dot 

We begin by considering a quantum system (e.g. quantum dot, molecule, etc) connected 
to two electrodes, left and right, with different chemical potentials. Each electrode is 
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Figure 1. Schematic illustration of quantum dot embedding. The electrodes are 
divided into macroscopic "environment" and buffer zone. The projection of the 
environment results into the Lindblad kinetic equation for the reduced density matrix 
of the buffer and quantum dot. Each buffer zone contains a finite number of discrete 
single-particle levels 

partitioned into two parts (Fig [I]): the macroscopically large lead (environment) and 
the finite buffer zone between the system and the environment. So the Hamiltonian of 
the whole system is 

U = H s + Hsb + Hb + H be + H E . (1) 

We assume that the environment and the buffer zones are described by the 
noninteracting Hamiltonians 

He = £ka a ka a ka, H B = £ba&\ a a ba- (2) 

ka ba 

Here S\. a denote the continuum single-particle spectra of the left [a = L) and right 
(a = R) lead states, a ka (ak a ) create (annihilate) electron in the lead state ka. 
The buffer zones have discrete energy spectrum E\, a with corresponding creation and 
annihilation operators a\ a and a\, a . The system Hamiltonian is taken in the most general 
form: 

H s = ^e s a\a s + H' s , (3) 

s 

where a\ (a s ) create (annihilate) electron in the single-particle state e s in the dot and H' s 
contains two-particle electron-electron correlations, and/or electron- vibration coupling. 
The buffer-environment and quantum dot-buffer couplings have the standard tunneling 
form: 

Hbe = ^2(v bka al a a ka + h.c), (4) 

bka 

H S b = 22{t s baal a a s + h.c). (5) 

sba 

Now we introduce an embedded system which consist of the quantum system itself 
and the buffer zones. We have recently demonstrated that if we take the buffer zones 
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sufficiently large the density matrix of the embedded system obeys the kinetic equation 
of Lindblad type. The technical details of the derivations and underlying approximations 
can be found in Appendix of [23] . Here we give only the sketch of the derivation with 
the emphasis on important physics relevant to our subsequent discussion. 

The starting point is the Liouville equation for the total density matrix x{t) in the 
interaction picture 

Xi(t) = -i[vj(t), X i(t)]- (6) 

Here the buffer-environment coupling Hbe is treated as an interaction Hamiltonian, i.e., 
T-L = h + Hbe and vi(t) = e tht H BE e~ lht . To derive the Lindblad master for the reduced 
density matrix of the embedded system, piit) = Tr EXi(t)i we take the trace over the 
environment in Eq. and make the following approximations: 

(i) The total density matrix can be factorized as Xi{t) — Piif)PE , where pE is density 
matrix of the environment taken in the equilibrium grand canonical ensemble form 
(Born approximation); 

(ii) The environment relaxation time is very fast, so we can use local-time (Markov) 
approximation for the reduced density matrix; 

(iii) The single particle states in the buffer zone propagate as free states 

e iht a ba e- iht = e'^a^ + 0(1/ VN) (7) 

where N is the number of discrete single particle levels of the buffer zone; 

(iv) Rapidly oscillating terms proportional to exp[z(£& a — ey a )] for Eb a ^ Sfa are 
neglected (rotating wave approximation). 

Under these approximations, the Liouville equation ([6]) reduces to a master equation 
for the reduced density matrix in Lindblad form. In the Schrodinger representation it 
can be written as 

^- = -i[H,p(t))+Up(t). (8) 

Here the Hamiltonian H includes the Lamb shift of the single-particle levels of the buffer 
zones 

H = H s + H SB + H B + J2 ^ba4 a a ba , (9) 

ba 

and the non-Hermitian dissipator is given by the standard Lindblad form 

npW = EE ( 2 W(*)4*„ - {L\ a ,L ha „ P {t)}). (io) 

ba A*=l,2 

The operators L haX and L ba2 are referred to as the Lindblad operators, which represent 
the buffer-environment interaction. They have the following form: 



J bal 



V^baldbcn L ba2 — \fTba2o\a- (H) 
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with T bal = 7 6a (l - f ba ), T ba2 = 'Jbafba- Here f ba = [1 + e^ Eb " ^] 1 and ^ ba (A ba ) is 
the imaginary (real) part of the environment self energy \v b k a \ 2 / (e ba — e^ a + i0 + ). 

The Lindblad master equation describes the time evolution of the open embedded 
quantum system preserving the probability and the positivity of the density matrix. 
Open boundary conditions are taken into account by the non-Hermitian dissipative part 
of Eq.®, Ilp(t), which represents the influence of environment on the buffer zone. The 
applied bias potential enters into Eq.® via fermionic occupation numbers f ba which 
depend on the temperature (/3 = 1/T) and the chemical potential fi a in the left and 
right electrodes. 



2.2. Liouville-Fock space 

Let us convert the Lindblad master equation JH} to a non-Hermitian field theory 
suitable for perturbative many-body calculations. To this aim we need to introduce 
the concept of creation and annihilation superoperators acting on the Liouville-Fock 
space [25j I2S1 [27J [22]. Our introduction of the Liouville-Fock space closely follows 
Schmutz work [25]. It is general and not restricted to the particular choice of the kinetic 
equation. 

Let {|^)} be a complete orthonormal basis set in the Fock space J 7 

^ \n){n\ = I, (n\m) = S nm . (12) 



nj\n) 



It is formed by particle number eigenstates \n) — a, . . . o] n |0), such that a^aj\n) 
Here |0) is the vacuum state and a], a,- are creation and annihilation operators for single- 
particle state j. Without loss of generality we focus on fermions, so we assume that a] 
and aj satisfy the canonical anti-commutation relations. 

The set of linear operators {^4(a^, a)} acting on T form a linear vector space, which 
is called the Liouville-Fock space associated with J 7 . We denote an element of the 
Liouville-Fock space by \A). The scalar product of two elements of the Liouville-Fock 
space is defined as 

(A 1 \A 2 ) = Tr(A\A 2 ). (13) 

In the Liouville-Fock space we introduce a complete orthonormal basis { |m, n) = 
\\m)(n\)}, which satisfies 

(mn\m'n) = S mm '5 nn ', \mn)(mn\ =1. (14) 

mn 

Here (mn\ = \mn)^ = ([\m)(n\] Ji \ = (\n)(m\ \ , and I is the identity operator in the 
Liouville-Fock space. Then, for an arbitrary element of the Liouville-Fock space we 
have 

\A) = J2 A mn\mn), (15) 
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where A mn = (m\A\n) = (mn\A). In particular, the identity operator I in Eq. (1121) 
corresponds to 

|I) = X>> n >- ( 16 ) 

n 

The scalar product of a vector \A) with (I\ is equivalent to the trace operation in the 
Fock space, 

(I\A)=Tc(A), (17) 

and for the density matrix we have (I\p) = 1. 

As was suggested by Schmutz [25] we introduce superoperators a, a through their 
action on the basis vectors \mn) 

dj \mn) — \aj\m)(n\), a,j \mn) = i{— 1) M ||m)(n|aj), (18) 

where \x = Ylj( m j + n j) = 171 + n - By analyzing the Hermitian conjugate of the matrix 
elements of a, a, we find 

dj \mn) = |aj|m)(n|), a] \mn) = i{— 1) M ||m)(n|a 3 -). (19) 

It follows from (|18[) and fTIT?]) that superoperators d, d^ simulate the action of a and 
on |m)(n| from the left, while a, simulate the action of a) and a on |m)(n| from the 
right. Here we would like to emphasize that our definition of tilde superoperators a, 'd) 
differs from Schmutz's definition by phase factors — i and +i, respectively. The reason 
for introducing these factors is that the so-called tilde-substitution rule (see bellow) 
becomes simpler. We also note that the alternative definition for superoperators is used 
in [27], where the "right" creation and annihilation superoperators are not Hermitian 
conjugate to each other. 

As follows from ffl8|) and ffT9l) . the superoperators ] obey the fermionic 

anti-commutation relations: 

{di,dj} = {d u a)} = Sij, (20) 

while other anti-commutators vanish 

{di,dj} = {ai,aj} = {di,%} = {d^a]} = 0. (21) 

It also follows from (USD and (JTSJ) that d |00) = a|00> = and the Liouville-Fock 
space basis vectors are generated from the vacuum 1 00) by application of the creation 
superoperators 

\mn) = {-ifdl i ...dij h ...ai\m). (22) 
Moreover, basis vectors \mn) are "superfermion" number eigenstates 

djdj \mn) = mj \mn), a^dj \mn) = rij \mn). (23) 

Using the definition of superoperators we can rewrite the identity ffTB]) in the 
following form 

|/)=expH^atat)|00). (24) 

3 
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Note, that because of the different definition of tilde superoperators, the obtained 
expression for \I) differs from Schmutz's analogous expression [25] by the phase factor 
(— i) in the exponent. From ( fl8|T9|) and ( 1241) we find that the superoperators a ] and a,j 
are connected to their tilde conjugate aj and aj by the relations 

a d \I) = -ioj a)\I) = -ia 3 \I). (25) 

For an operator A = A(a', a) given by the power series of creation and annihilation 
operators we define two superoperators 

A = A(a\a), A = A*(a\a). (26) 

Here, the * means the complex conjugate of the c-number coefficients. The relation 
between non-tilde and tilde superoperators is given by the following tilde conjugation 
rules 

(c 1 A 1 + c 2 A 2 y=clA 1 + c* 2 A 2 , (A 1 A 2 y=A 1 A 2 , {AJ= A. (27) 

Applying tilde conjugation to |mn) we find 

\mnY= (+i) tl2 Iran), (28) 

where p = m + n. Therefore |IJ~ = \I), i.e., \I) is tilde-invariant. Generally, if 
A = A(a^,a) is a Hermitian bosonic operator then |AJ~= \A). 

According to the definition of the superoperator A, if A = J2mn J ^-mn\ m )( n \ then 
A = ^2 mnk A mn \mk)(nk\ and we obtain 

\A)=A\I), \AT=A\I), (29) 

\A 1 A 2 ) = A 1 A 2 \I) = A 1 \A 2 ). (30) 

Therefore, the expectation value of an operator A = A(a\a) in the state with the 
density matrix p = p(a\ a) can be calculated as the matrix element of the corresponding 
superoperator A = A{a\a) sandwiched between (I\ and \p) = p\I) 

(A) = Tr(Ap) = (I\Ap) = (I\A\p). (31) 

Using fl25|) we can show that the following tilde-substitution rule is valid 

A\I)=a A fi\I). (32) 

Here a a = +1 if A is a bosonic operator and a a = —i if A is a fermionic operator. 
Moreover, taking into account that non-tilde and tilde fermion superoperators anti- 
commute we find that 

A 1 \A 2 )=iA\\A 1 ) ) (33) 
if both A\ and A 2 are fermionic operators, and 

ia \A 2 ) = <jaA \Ai) (34) 
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otherwise. It should be noted that Schmutz tilde-substitution rule |25J is cumbersome 
and it takes the simple form like ( 132]) only if all terms in the power series of A(a\a) 
have the common quantity m — n. Here m(n) is the number of creation (annihilation) 
operators. 

The general prescription to obtain equation for \p(t)) from the kinetic equation 
for p(t) is the following. First, we transform the kinetic equation for p = p(a),a) 
into the kinetic equation for p = p(a',a) by formally replacing all operators a', a by 
superoperators a', a. Then, we multiply the kinetic equation from the right on vector 
| J) and use (|32|) - flM|) to convert the kinetic equation to the Schroedinger-like equation 
for the vector \p(t)) = p{t) 

Ip(*)> = L(a^,a,a i ,a) \p{t)), (35) 

at 

where L is the Liouvillian which depends on both non-tilde and tilde superoperators. 
Particularly, the Liouvillian for the Lindblad master equation (j8]) becomes 

L = H - H - ij^^ba, (36) 

ba 

where 

riba =(r 6a i — T ba2 )(al a a ba + a\ a a ba ) 

— 2i(T ba ia ba a ba + r& a2 a 6Q ,a bc J + 2T ba2 . (37) 

In derivation of (1361) and (1371) we took into account that p = p(a\a) is a bosonic 
superoperator which commutes with all tilde superoperators. Due to the Lindblad 
dissipators, the Liouville superoperator ( |36l) is non- Hermit ian. In addition, as \p) is 
tilde- invariant, the Liouvillian obeys the property [Ly= —L. 

Taking the time derivative of (I\p(t)) = 1 we find that (I\ L = 0, i.e., (I\ the left 
zero-eigenvalue eigenstate of the Liouvillian superoperator. Since also (J| is the vacuum 
for a] — ia,j and a] +ia,j superoperators, it is appropriate to call (J| left vacuum vector. 
For the electron transport problem we focus on nonequilibrium steady-state where the 
current through the quantum dot is given by 

(J a ) = Tr(J a poo) = (I\ Jalpoo). (38) 

Here J a is the current superoperator, and the stationary, steady-state solution of ( J35l) . 
Ipoc), is the right zero-eigenvalue eigenstate (right vacuum vector) of the Liouville 
superoperator 

£|Poo> = 0. (39) 

In the next section, we show how one can find Ipoo) perturbatively starting from the 
free-field approximation for nonequilibrium density matrix. 
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3. Perturbative calculations of the steady state density matrix and electron 
current 

3.1. Nonequilibrium many-body perturbation theory 

Let us make the important remark on the notation use in the rest of the paper: only 
creation/annihilation operators written with letters a,d (such as for example and 
a) ba ) are related to each other by the Hermitian conjugation; all other creation eft ,b\rf 
and annihilation c, b, 7 operators are "canonically conjugated" to each other, i.e., for 
example, does not mean (c)' although cc 1 ± (ftc = 1 (± - bosons/fermions). We will 
also use the same notation for the non-tilde superoperators a\ and a,- as the ordinary 
operators aj and aj bearing in mind that all operators acting in the Liouville-Fock space 
are are superoperators. 

We start by rewriting the Liouvillian f )36|) as 

L = L (0) + U, (40) 
where Z/°) is the quadratic unperturbed part of L, and 

U = H' s - H' 8 (41) 
is a perturbation. Then using the equation of motion method 

[4l (0) ] = 

[c n ,L^] = Q n c n , (42) 

we exactly diagonalize|22j in terms of nonequilibrium quasiparticle creation and 
annihilation operators: 

L<°> =J2^nclc n -Q*Jicn). (43) 

n 

Here c^ n ,c an are obtained from eft an , c an by the tilde conjugation rules. 

The nonequilibrium quasiparticle creation and annihilation operators are connected 
to aft , a, 'aft, a by canonical (but not unitary) transformations: 



s ba 

y^(V>n,A + i<Pn,sbl) + ^(tpufiabba + i¥n,ba b ba)> ( 44 ) 



where 



(k r , . . .. 

s ba 



b\ = a\- ia s , b, = a.,, bL 



's "'s """»> "8 "ha ~ a ba ^ a bo 

bba = (1 — fba) a ba + ifba&\a' 



Nonequilibrium quasiparticle creation and annihilation operators obey the fermionic 
anti-commutation relations. In particular, from {c n ,cL,} = b~ nn i and {c„,c n /} = we 
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find the following orthonormality conditions for amplitudes 

^ i ! n,s' l Pn',s + ^ i'nM'^ri 1 ,ba = &nri' , 
s ba 

s 

+ ^(^nMVl'M - VnM^l'.ba) = °- ( 45 ) 
ba 

By construction, (I| is the left vacuum for c^, operators. The vacuum 
for c n , c n operators, \poo), is automatically the zero-eigenvalue eigenstate of the 
unperturbed Liouvillian L^°\ i.e., it is the steady state density matrix in the zeroth-order 
approximation: 

i (P) |p£ ) > = 0, (J|pS> = l. (46) 

In other words, the the zeroth-order density matrix is the density matrix which does 
not contain nonequilibrium quasiparticle excitations. 

Now we introduce the continuous real parameter A, which will be set to unity in 
the end of the calculations, 

L = L (0) + XL' (47) 
and expand the exact steady state density matrix in powers of A, 

|Poo> = £A^)>. (48) 

Substituting (148]) into Eq. (139|) . we obtain equation for the pth-order correction to the 
zeroth-order density matrix: 

L \pto) = -L'\p£-% (49) 

or |/?So ) = (— Lq l L') p \p ( £}). Here, V is expressed in terms of the nonequilibrium 
quasiparticles. Thus, starting from \p^) we can find any pth-order corrections to 
it. In addition, (/| p&) = for p > 1 since \pa$) contains excited nonequilibrium 
quasiparticles. 

To calculate the current through the quantum dot we express the current 
superoperator 

J a = -i t sba (al a a s - a\a ba ) (50) 

bs 

in terms of nonequilibrium quasiparticle creation and annihilation operators and 
compute its expectation value with respect to (I\ and \poo)- As a result we get the 
following expansion 

j a = Y. XPJ « ] - ( 51 ) 

p=0 
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Here, Ja^ is zeroth-order current for the system without interaction 

J^ 0) = -21m t sba ip nM Lp njS , (52) 

bsn 

and J« p ^ is the pth-order correction to it 

J<W = -21m ^ ta^l^.F^l (53) 

where _Fm2 is the expansion coefficient in 

|pS>=iE F ™^|p£ } > + ..- (54) 

mn 

and F„f n = (i*nm)* as follows from ) = |poo J~- Thus, the problem of computing the 
pth-order correction to the unperturbed current is reduced to finding F mn by solving 
Eq. (SHD- 

Using the same method we can obtain perturbative expansion for the population 
of a quantum dot single-particle level 

n s = {I\a\a s \ Poo ) = Y,^\ (55) 

where 

n mn 

The anti-commutation condition {6 S ,& S } = imposes the constraint on the amplitudes 
from which follows that is a real number. 



3.2. Electron-vibronic coupling 

As the first application of the method we consider the Hamiltonian Hs which describes 
one electronic single-particle level coupled linearly to a vibration mode (phonon) of 
frequency u (so-called local Holstein model) 

H s = e^a^a + uj Q d)d + K(va{jv + d). (57) 

For simplicity we assume that the tunneling matrix element in Eq. (jSJ) is real number 
t independent of indices a and b. The electron spin does not play any role here, so 
we suppress the spin index in the equations in this section. Replacing k by \k, we 
arrive to perturbation expansion of the steady state density matrix \poo) with respect 
to electron-vibronic coupling 

| Poo ) = EA p |p^). (58) 

p=0 
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To find the zeroth-order density matrix \poo), we diagonalize the fermionic part 
of L/(°\ The resulting creation and annihilation operators have the form (144p . and 
amplitudes if), <p satisfy the following system of equations 

£0^n ~tJ2 ^n,ba = ^n^n 

ba (59) 

Ebatpn.ba tlf) n ^n^Pn,bat 

(Eq — Cl n )(f n — t X] fn,ba = t X, fbalf>n,ba 

ba ba (60) 

(Eba — ^n)fn,ba ~ tip n = —tfb a lf> n , 

where Eb a = Sb a — ijba- The solution of eigenvalue problem ( 159]) yields the spectrum 
of nonequilibrium quasiparticles, Q n and — f2*, as well as if) amplitudes which should 
be normalized according to Eq. f )45|) . To find (p amplitudes we must solve linear 
equations ( l60l) . 

Furthermore, let be the number of vibrational quanta with frequency Uq at 
temperature 1/(3, i.e., N u = (exp(/3co> ) — l) -1 . When k = the density matrix is 
factorizedas |p£ } ) = |^> / |pg?) 6j 

(/| dtd |pg?> = (61) 

It is convenient to introduce new phonon operators 

1 =(l+N u )d-N u $, 

^ = S -d (62) 

and their tilde conjugated partners, such that (I\ 7^ = (I\ 7^ = and 7 \p^} = 7 \poo) = 
0. Then the quadratic part of the Liouvillian is diagonal in terms of introduced operators 

L<°> = ^(O n 4c n - Q*^ c n ) + ^(^7 - 7 f 7), (63) 

and the vacuum for c n , c n , 7, and 7 operators is the the zeroth-order approximation 
for the density matrix, \p$ }■ For the unperturbed current we have 

J^ = -2tImJ2^n, (64) 

bn 

while pth-order correction is 

= -2tImJ2r nM ^ m F!Zl (65) 

bmn 

(v) 

To find Fmn, we rewrite the perturbative part of Liouvillian in terms of operators 
c n , 7, etc.: 

U = J2{ Kb f + & + Lffi(7 + 7)Rc„ - t.c.} 

~~ * ^ ] [-^mn7^ ~~ (^iii) 7^ "I" ^rnni'l 7)] C m,^n 
mn 

-< E L ™(7 f - ^) c ^« + ^ (0) (7 f - 7 f ), (66) 
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where coefficients LW are 



L mn = K[{lp m ~ (f m ) + N^ m ]^ n , 

L$ n = n[ip m + N w 4> m ]4> n , L$ n = K^m^n, 

L mn = «[(^m ~ <Pm)<Pn + N u (lj; m (p* ~ ¥W0»)] 

L mn = K bPmfn ~ fm^n] , L mn = VVnVC (67) 



and 

,(o) 



a' =(I\a ] a\p^)=^2ip n ip n (6£ 



is an unperturbed electron level population. The notation 't.c' in equation (1661) means 
the tilde conjugation (i.e., c m — > c m , 7^ — >■ 7, Lmn — >■ (Lmn)*, etc.). Then, substituting 
Eqs. ( |63l 1661) into ( 1491) we obtain the following general expression for Fmn 

Fmn = ~ q _ O* I [^kn ^ + i^nk ) ] 

-E^SVK^V + ^l -2L$ n W^}, (69) 

fc 

where Z™ and are coefficients in the expansion 

|pg> = {w«( 7 t + ^t) 

+ < + (^S^R^ + • • •} Ip£ } >- (70) 

mn 

Thus, to find pth-order correction to the current we need first compute Zmn and 
W(p-i). This can be down using the same method as used to find Fmn- As a result, Zmi 
and are nonvanishing only for odd p. Therefore, only even powers of p contribute 
to the current expansion as it should be for the considered model. It is interesting to 
note that the term W^irf + 7^) \p<x) is associated to the momentum transfer from the 
electronic current to the quantum dot vibrational mode (current induced translational 
motion of the dot) whereas Zmn^c^c^ and (Z n p m )*^ c^ch ) correspond to the 
current induced heating and cooling processes respectively. 

As an example we give here explicit expressions for the first order perturbation 
theory W {1) and Z { „ 1) ■ 



J mn- 



W^ = - — } ZW=—-^ . (71) 

_____ /fj\ /fj\ 

Combining Eqs. ( 1711) and (1691) . we find F mn . Then inserting F mn into ( 1651) we derive 
the second-order perturbation theory correction to Ja ■ This correction consists of 
two parts: the first part is proportional to n^ -*, so it is the Hartree term, while the 
remaining part is the Fock term. In section 13.41 we also verify these definitions by 
comparing Hartree and Fock terms obtained within the presented approach and the 
exact ones given by NEGF formalism. 
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3.3. Electronic correlations 

As a next example we consider electron transport through one spin-degenerate level 
with local Coulomb interaction 

H s = £ o^2 n ° + Un^. (72) 

(7 

Here n a = a\a a is the number operator for electrons with spin a in the quantum dot. 
In what follows, we again assume the tunneling matrix element is independent of a, b 
as well as spin a, i.e., 

H S B = -tJ2( a L« a «+ h - c )- ( 73 ) 

aba 

We also assume that energy levels in the leads are spin-degenerate. 

Since the quadratic part of the corresponding Liouvillian describes electron 
transport through noninteractiong spin-up and spin-down levels, it is diagonalized by 
the same method as in the previous example. As a result we obtain 

L<°> = ^(O n 4 n C CTn - K^nCan). (74) 
an 

The vacuum of c an and c an operators, \p6o), is the density matrix in the zeroth-order 
perturbation theory and 

= -4tIm^Vw,a^n (75) 

lm 

is the corresponding current. 

To find pth-order correction to fl75l) . 



J« = -4ffm Kba^mF^, (76) 



bmn 



we rewrite V = U(n^n^ — n^n^) in terms of nonequilibrium quasiparticles: 

u = X){tf2 } (4*cai - t.c.) + 



akl 



fc/mn C l t C Zj. C »H C "t 1 ' * ' 

+ L (2) c f c f # # 

~ C k 1 - C l- t C m l c n i — c i i c l i c m 1 .C„ t ) 
+ * [-^Lmn( C l t C i\^mj. C nt + ^k^l^m^ C n^) + t.C.] 

+ « [ L £L( c I t ^ c ^ c "t + cl^t^t^J + Lc -] }• ( 7T ) 
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Here and K^ are given by 

(0) -</|4o ff |p£)> = X;^», (78) 



n. 



while the coefficients L^j mn are listed in 



for Fmn 



Now, substituting Eqs. (I74|77l) into Eq. fH9|) we find the following general expression 
= Ik W 5 , 

mn q I 11 mn u pf 

+ V[#F M -f/i' (1) F (H Vl _ VV( 3 ) F {p ~ l) 

"i - / >[ ini in v to im / J / , mnij ji 

i ij 

'^^[^rnijk^k'jn? ~ (^nijk^kjmi)] }' (79) 
(P) 

mi uena anu kjt 

roo 



ijfc 

where 5 p i is the Kronecker delta and G^ mn is a coefficient in the expansion 

» 



} ) = {E ^4*4*4^ + • • •} IpS>- (so) 



In turn, the equation like (1791) can be derived for Cr 



klmn 

klmn' 



The exact first-order perturbation theory correction to |p6o ) is 



mn^ffm^ffn 



, ^(1) „t „t^t \| n (0)\ /oi\ 

^fc«mn c ffc Ji 1™ I™/ I Poo /> l 01 J 



amn 



where 



klmn 



K m , N t (2) 



Ei(l) _ -"-mn ynr(l) ^klmn 

fim-n*' klmn n k + ^-n* m -n* n 1 J 

Inserting Fmn into (176]) we get the first-order perturbation theory correction Ja to the 



current (I75|) . This correction is proportional to n^ ) and below we will show that it 
corresponds to the first-order Hartree term obtained with NEGF formalism. 

Here we note, that in [23] we applied perturbation theory to the Anderson model 
starting from the nonequilibrium Hartree- Fock approximation, i.e., V was normal 
ordered and did not contain quadratic terms. Therefore, in [23] the mixture of two 
quasiparticle configurations to |pco) vanished and the first-order perturbation theory 
correction to the current was zero. 

To find the second-order correction to Ja we insert ( 1B2I) into ( 1791) . This yields 

F& = \ {y \K {1) F {1) - (K {1) F {1) r] 

mn q _ I / j I mi in V to im ! J 

Ei- (3) p( 1 )_^rr(5) _,-r(5) r (i) /oo\ 

mnij ji / j L miik kini \ nijk kjmil J J" v od / 

ij ijA; 
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Figure 2. The second-order perturbation theory correction to the current for the local 
Holstein model: Hartree term. 

( < ?\ 

Now, with the help of the obtained expression for Fmn and Eq. ( 1761) we get the second- 
order perturbation theory correction to Ja . 

3.4- Comparison with Keldysh NEGF perturbation theory 

Let us now compare the results obtained with the present approach with those calculated 
with the help of Keldysh NEGF. For the case when the coupling to the left lead is 
proportional to the coupling to the right lead, the electron current through the quantum 
dot can be computed directly from the retarded dot Green's function, G r (ui), using the 
well known Meir-Wingreen formula [28]. For the considered models, assuming that the 
left and right leads are identical, Tl ^(uj) = 0.5T(u>), this formula takes the form 

J= h\ duJ ^ L ^ ~ fR(u)]T(uj)ImG r (uj). (84) 

Here s is the spin degeneracy of the considered models: s — 1 for the model with 
electron-vibration coupling and s = 2 for the model with electron-electron interaction. 
We will use the wide-band approximation for the electrode, so the imaginary part of the 
electrode self-energy, which is responsible for level broadening, is energy independent, 
r(o;) = T, while its real part vanishes. 
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Figure 3. The second-order perturbation theory correction to the current for the 
local Holstein model: Fock term. 



The retarded Green's function is the solution of the Dyson equation 
G r (u) = G r (cu) + G r (u)i: r (u)G r (u), 



J5) 



where Gq(oj) = (u — Eq + iT)~ l is the noninteracting retarded Green's function and 
E r (o;) is retarded self-energy evaluated in the presence of electron-electron or electron- 
vibration interaction. Expanding S r (w) with respect to electron-electron or electron- 
vibration coupling, X r (w) = -^^pt^); we obtain perturbative expansion of G r {oS) 



J3=l 



and consequently of the current 



X 



p=l p=0 



(86) 



Here is the current through the system without interaction given by the standard 
Landauer formula. 

In [22] we have shown that for the current through a system without interaction, 
the exact agreement between NEGF and kinetic equation approach can be achieved 
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Figure 4. The first-order perturbation theory correction to the current for the 
Anderson model. 



by increasing the density of states in the buffer zones. Below we show that this is also 
true for correlated electronic systems. 

In what follows, in the calculations based on the kinetic equation we will assume 
that iV single-particle levels in each buffer zone are evenly distributed within the energy 
bandwidth [E m i n ; £ max ] = [—10,10]. This bandwidth is much larger than any energy 
parameter in the system, so it corresponds to the wide-band approximation used in 
NEGF calculations. The tunneling coupling strength t is computed from T = 2nr)t 2 , 
where rj = N/(E max — E m[n ) is density of states in the buffer zone. We note here that the 
main approximation in the derivation of the Lindblad master equation (jH|) is that the 
single particle states in the buffer zone propagate in time as free states (j7|). It is evident 
from Eq.flZJ) that the larger the buffer zone, i.e. the larger the density of states r], the 
better this approximation. This will be also explicitly demonstrated in the numerical 
calculations below. The parameter 7 in the Lindblad operators is chosen to be 7 = 2Ae, 
where e is the energy spacing between the energy levels in the buffer zone. In addition, 
although it is not necessary, we use a symmetric applied voltage, [Il,r = ±0.51/, and 
the temperature of the electrodes is zero, T = 0. 

At the beginning, we consider the system with electron-vibronic interaction and 
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Figure 5. The current through the Anderson model computed by taking into account 
the first- and second-order corrections. 

compare the second-order correction to the current obtained in the section [372] with that 
calculated using NEGF formalism ( )86l) . We use the following model parameters of the 
Hamiltonian fl57|) : k = 1.0, Uo = 1.0. 

In NEGF formalism the second-order correction to the current arises from the 
retarded self-energy ££ which contains contributions from Hartree and Fock diagrams, 
= £# + Ti r F . The Hartree self-energy is [8] 

TT E ( U ) = -— n<°>, (87) 

where is the electron level population in the zero-order approximation 



n 



(o) = il + 



2% J (uj-soY + T^ 



The expression for the Fock self-energy is more complicated and can be found elsewhere 
(see, for example, [29]). 

In Figs. [2] and [3] we compare Hartree and Fock second-order corrections to the 
current obtained within our approach with different size N of buffer zone and the exact 
ones. The corrections are shown as functions of the level energy, Eq, for two values of the 
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applied voltage V and broadening T. It is evident from the figures that the difference 
between exact and Lindblad equation based results become negligible as we increase the 
leads density of states in the buffer zone. The reason is that increasing the number of 
single-particle state in the buffer zones we make the approximation (iii), under which 
Lindblad master equation (IE]) was derived, more justified. The deviation of the results 
obtained from the Lindblad kinetic equation and NEGF becomes smaller at the larger 
applied voltage or T. 

Now we compare first-order corrections to the current for the Anderson model. We 
put U = 1.0 for the strength of the Coulomb interaction. Within the NEGF formalism 
the first order correction is solely due to Hartree diagram and it is 



The results of numerical calculations are shown in Fig. H] for different values of F 
and applied voltage V. As we can see the results of the Lindblad equation approach 
converge to the exact results with increasing value of N and the convergence is faster 
for larger values of applied voltage and T. 

In Fig. [5] we show the current through the Anderson model computed by means 
of Lindblad equation by taking into account the first- and second-order corrections. 
We take N = 1600, so the obtained results correspond to NEGF ones. As we can 
see from the figure, the first- and second-order contributions shift the maximum of the 
current towards the symmetric point 6q = —0.5U. The first-order correction increase 
the maximum current, while the second-order correction acts in opposite direction. We 
also see from Fig. |5] that for a given U the relative value of the first- and second-order 
corrections show little dependence on the applied voltage V. In contrast, in [23] we have 
observed that nonequilibrium post-Hartree-Fock electronic correlations play important 
role at larger applied voltages and, as a result, the second-order correction to the 
current become more pronounced with increasing V. This is due to the difference in the 
structure and spectrum of nonequilibrium quasiparticles. The quasiparticle spectrum, 
both tp and if amplitudes depend on the voltage in the post-Hartree-Fock perturbation 
theory [23], whereas in the present work the voltage enters only into (p amplitudes of 
the nonequilibrium quasiparticles through Fermi-Dirac occupation numbers of the buffer 
states. 

4. Conclusions 

We developed nonequilibrium many-body perturbation theory for steady state density 
matrix and electric current through the region of interacting electrons. Our approach is 
based on the super-fermion representation of quantum kinetic equations. We considered 
an quantum dot connected to the reservoir through the buffer zone (so-called embedded 
quantum dot). The Lindblad type kinetic equation were obtained for the embedded 




(89) 



where the population 



is given by Eq. (J8"8~j) . 
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quantum dot and the kinetic equation was converted to the non-Hermitian field theory 
in Liouville-Fock space via the tilde conjugation rules. The free- field state was defined as 
vacuum for non equilibrium quasiparticles and this state describes the ballistic transport 
with the results equivalent to the Landauer formulae. We applied the nonequilibrium 
perturbation theory to compute corrections to nonequilibrium quasiparticle vacuum 
for the system with electron-phonon and electron-electron correlations. The exact 
agreement with the Keldysh NEGF perturbation theory was observed for inelastic 
electron current through quantum dot. 
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